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ABSTRACT 

We present analytical solutions for the integral distribution of arbitrary bursting or steady source counts 
as a function of peak photon count rate within Friedmann cosmological models. We discuss both the 
standard candle and truncated power-law luminosity function cases with a power-law density evolution. 
While the analysis is quite general, the specific example discussed here is that of a cosmological gamma- 
ray burst distribution. These solutions show quantitatively the degree of dependence of the counts on the 
density and luminosity function parameters, as well as the the weak dependence on the closure parameter 
and the maximum redshift. An approximate comparison with the publicly available Compton Gamma 
Ray Observatory data gives an estimate of the maximum source luminosity and an upper limit to the 
minimum luminosity. We discuss possible ways of further constraining the various parameters. 
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1. INTRODUCTION 

Gamma-ray burst (GRB) sources appear to be distributed isotropically down to the current limit 
of spatial resolution (Meegan et al. 1992, Fishman et al. 1994, Meegan et al. 1994), suggesting a 
cosmological origin. The brightness distribution, or the integral number counts of GRB as a function of 
peak photon count rates have been discussed in this context by, among others, Mao & Paczyhski (1992), 
Dermer (1992), Piran (1992) and Wasserman (1992). Statistical fits to this "log N - logC" distribution 
have been carried out recently by Wickramasinghe et al. (1993), Horack, Emslie & Meegan (1994), Cohen 
& Piran (1995), and Emslie & Horack (1994). Some of the questions addressed in these references include 
the degree of sensitivity of the GRB count rates to the details of the cosmological model, the density 
evolution and the luminosity function of the bursts. 

Except for the simpler standard candle case, these calculations have required extensive numerical 
integrations in order to explore the parameter space fully. These numerical calculations have shown 
that the fits are not very sensitive to the exact cosmological model, and that bounds on the density 
evolution and luminosity function are fairly general in nature. The role of each individual parameter in 
such numerical fits can only be inferred by varying them one at a time and recalculating the fit in each 
case. The analysis relies therefore on extensive sets of graphs and tables. A physical interpretation of 
the reason why the fits vary as they do as function of the different parameters would be greatly aided if 
one had analytic expressions for the quantities being fitted. These could, furthermore, greatly reduce the 
amount of numerical work involved, while making it easier to explore the effect of alternative distribution 
forms or evolution laws. 

In the present paper, we derive closed form analytical solutions for the integral burst number count 
rates for different cosmological models, including a source density evolution proportional to a power of 
the expansion scale and a power-law luminosity distribution with upper and lower limits. These analytic 
expressions provide direct quantitative measures of the relative contributions of these three different 



physical factors, and show explicitly the reasons for which the comparisons of the models to the data are 
largely insensitive to the specific values of the closure parameter Sl and the density evolution power law 
index. The analytic solutions show explicitly the functional similarities and differences between a standard 
candle interpretation and a power-law luminosity distribution case. For a luminosity distribution, the 
analysis provides a direct estimate not only of the luminosity function slope, but also of the maximum 
burst luminosity, as well as lower limits for both the ratio of maximum to minimum burst luminosities 
and the total burst rate per galaxy. In this paper we discuss approximate fits-by-eye of the analytical 
solutions to the observed log N — log C curves data, and the general behavior of the analytical solutions. 
These solutions are expected to be very useful in checking numerical integrations, and also in performing 
detailed statistical fits, which we intend to discuss elsewhere. 

The analysis considered here is fairly general and may apply to other phenomena as well, e.g. supernova 
outburst rates in external galaxies, extragalactic gravitational wave or neutrino burst events, radio galaxy 
flare-ups, etc. These are also phenomena where luminosity function and evolutionary effects can play a 
large role, and are not treated here. It can also apply to steady (non-bursting) sources, if one replaces 
the parameter D by D ss = D + 1, where D has the same meaning as in eq.(2). Our discussion here is 
developed for the specific case of a cosmological gamma-ray burst distribution. 

2. BURST RATE AND PHOTON COUNT MEASUREMENTS 

The number of generic transient events (henceforth bursts, or GRBs) observed per year out to some 
redshift z can be calculated if one knows the local burst rate density as a function of redshift. Since 
one is interested in observable bursts, their intrinsic luminosity must be taken into account, as well as 
the fraction of it which is actually observable by a particular instrument. Consider a GRB with a peak 
bolometric photon emissivity C (photons per second) at a redshift z. A K-correction factor will take into 
account the spectrum of the source and the fact that the detector measures only within a limited band of 
frequencies v\ and i>2 (Mao & Paczyhski 1992). If the photon spectrum is characterized by a power law 
C v ~ v a ~ 2 , where a is the index for the power-per-decade spectrum below (generally different below and 
above several hundred keV, c.f. Higdon & Lingenfelter 1990), then the effective measured "bolometric" 
peak emission rate will be 

£= C v dv={l + z) (a -^C K jC k = C v dv. (1) 

•J Ui(l-\-z) J V\ 

For a ~ 1 (typical of GRB spectra below about 0.5 MeV, where most of the measured photons are (e.g., 
Band et al. (1993)), the K-correction is not present. In what follows, except for the short note in §2 
below eq.(9) we will assume a = 1 and Ck = £■ 

The typical photon emission rate C and the number density of bursts may depend on the redshift z. 
We denote by $(£,z)cLC (in Mpc -3 yr _1 ) the number of GRB per year per unit proper volume with 
emissivity between C and (C + dC). This may be parametrized as 

$(£,z) = $(£)(l + z) c , 0<Z?<4, (2) 

where $(£) = $(£, 0), and (1 + z) D is used to define the density dependence. D = 3 would correspond 
to a constant comoving density, while for D > 3 (D < 3) the comoving density decreases (increases) with 
time. (The restriction on the values of D is a matter of convenience only, other values being in principle 
allowable too. In this paper, for simplicity, we restrict ourselves to integer values < D < 4). The total 
number of GRB per Mpc 3 per year at redshift z is then 

n(z)= [ ™" $(£)rf£ = (1 + z) D [ $(jC)djC = (l + z) D n , (3) 



where n = n(0). We assume that < C m i n < C < C max < oo. 

The volume integrations are most conveniently carried out in terms of the comoving radial distance x 
or the dimensionless conformal time rj. The relationship between these quantities and the redshift is (e.g. 
Weinberg 1972) 

X = Vo-V, (l + z) = R /R( V ), (4) 

where rj and R = R(r) ) are the present conformal time and present expansion scale factor. The 
total number of sources located at comoving distances between x an d X + i s ( e -g- Meszaros & 
Meszaros 1988) ^n{x)R?{r\)f 1 {x)dx ) where f(x) is either sinhx, X; or sm X f° r an open, flat or closed 
universe. Due to time dilation, the number per year of bursts at the observer position x = is (1 + z)~ x 
times smaller (Mao & Paczyhski 1992), or (1 + z)~ 1 4iTn(x)R 3 ('ri)f 2 (x)dx- Hence the total number of 
GRB observed per year at x = out to comoving distances < xi (the integrated GRB count rate) is 

N'(< X i) = ^n Rl r\R(il)/Ro) A - D f(x)dx- (5) 



In the Euclidean case, the bolometric peak count rate C (in photons cm 2 s 1 ) observed from a GRB 
at distance d with bolometric peak photon emission rate C photons per second is given by C = C/{Ai:d 2 ). 
The corresponding expression for the Friedmann model case for a GRB at comoving distance x is (aside 
from possible K-corrections) 

r = — = C = CR(r]) ((\) 

47Trf2 47Trf2(l + z) AttR^X 2 ' ( ' 

where d p = R X is the proper distance (Weinberg 1972) out to the redshift z corresponding to the 
comoving distance x, (1 + z ) is the redshift factor accounting for cosmological time dilation, and we have 
defined a photon "count distance" d c = d p (l + z) 1 ! 2 related to the more usual "luminosity distance" 
d, = d p (l + z) by d c = d,(l + z)- 1 ' 2 . 



3. STANDARD CANDLE DISTRIBUTION IN A FLAT UNIVERSE 



The simplest luminosity distribution is one where GRBs are standard candles, $(£) = n 8(C — C ), 

where J £ 8(C — C )dC = 1. For the Einstein-de Sitter flat universe, the function fix) = X; an d one 

can take r\ = 1; R = 2c/H ; R(r]) = R r\ 2 = R (l — x) 2 , where c is the velocity of light and H is the 
present Hubble parameter, while the count distance is 

d c = d p (l + z) 1 ' 2 = d,(l + z)- 1 ' 2 = R oX (R /R(ti)) 1/2 = R (VlT~z - 1) . (7) 

In this special (flat universe) case, the counts received from x are therefore C = (£ /(4:'irR 2 ))((l — x) 2 /x 2 ), 
and the number of bursts observed per year originating from comoving distances up to xi is 

N(< xi) = ^n R 3 j X \l - x) 2(4 - D) X 2 dx, (8) 
Jo 

assuming that the burst spatial distribution sample is brightness limited (or rather peak count limited, 
Xi *C Xmax', for a redshift limited sample, see Appendix B). Considering integer values < D < 4 and 
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substituting xi = [1 + Roi^C/Co) 1 !' 2 ]- 1 , eq.(8) gives a closed expression for the integrated burst count 
rates per year 

3/2 (8-2D) 



N(> C) = ^n -^— 2 (1 + (C/VkCRI)) 1 '*)- 3 [J2 a * C 1 + (^Ci? 2 /^) 1 / 2 )"^ 

^ ' k = 



* ~ (Jfe + 3) k\(% - 2D - k)V { ' 

Equ.(9) defines the function N(> C) of the independent variable C with three parameters n , C , and 
D. (This expression and the definition of a k is valid also for half-integer values of D = 1/2, 3/2, 5/2, 7/2, 
besides the integer values D = 0, 1, 2, 3, 4 used here as examples). For non-integer D the integral in equ.(8) 
can also be calculated analytically in particular cases, and numerically in general. The analytic expression 
(9) with integer values of < D < 4 delimits the range of behavior expected for the integrated burst 
count rates over most cases of interest. For the special case D = 3 (constant comoving GRB density), 
equ.(9) reduces to a form previously derived by Mao & Paczyhski (1992) and Wasserman (1992). Note 
that, if a ^ 1, D in equations (5) and (8) can just be replaced by D e jj = (D + a— 1), and in equ.(8) the 
integration can be done very simply for integer and half-integer D e jj. In this case, one needs in addition 
to use C = (1 + z)"~ 1 Ck in eq.(6), and hence here the relation between \ an d £-k = £o is given for 
any C by C = (£ /(47ri? )) 2 ((l — x) 4_2 "/x 2 )- Calculating \ from this relation, and substituting it into 
the integrated version of equ.(8), one obtains the analogue of equation (9). In what follows, however, we 
continue to use for simplicity the case a = 1. 

In the limit of very large and very small C the integrated burst counts tend to 
lim N(> C) = (47r/3)n £o /2 (47rC7)- 3 / 2 ; 

C— >oo 

lim N(>C)= (47r/3)n i^ [3/(11 - 2D) -6/(10 - 2D) + 3/(9 - 2D)]. (10) 

(The second limit can be obtained from equ.(8) by taking for the upper limit xi = 1 corresponding 
to z ^ oo; the square bracket in (10) is identical to the Ad defined in eq.(21)). This means that 
the theoretical log#(> C) vs. logC curve is a monotonically decreasing function, with d(logN(> 
C))/d(logC) monotonically decreasing from to — | as C runs from to oo (as log C runs from — oo to 

oo). In reality, of course, N(> C) as a function of C is defined only for C > Cth > 0, where Cth is the 
observational threshold peak photon count rate. 

The relation (9) can be used for a limited test of the cosmological interpretation of GRB, under the 
assumption of a flat universe in the brightness limited standard candle case. Since %i = 1 — (1 + z) -1 ! 2 , 
we can rewrite equ.(9) in the form 

a r 3 / 2 

N{> C) = —J^cj^ qiz) = 

a r 3 l' 2 {S ~' 2D) 

[(! + z )- 3 ' 2 E (!"(! + ^)" 1/2 )1 > (") 

where < q(z) < 1; q(0) = 1 and is defined in eq.(9). The function q(z) (the square bracket in 
equ. 11) quantifies the degree of departure of the curve N(> C) vs. C from the Euclidean homogeneous 
unbounded case, for which N(> C) ~ C~ 3 I 2 . 
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It is instructive to carry out a rough visual comparison of the data against the simple theoretical 
model considered in this section. For any C > Cth one has from the BATSE data an observed value of 
N(> C), whereas from the extrapolated -3/2 relation one obtains for any C > Cth a second "Euclidean" 
value Ne(> C). The ratio N(> C)/Ne(> C) must be given, under the standard candle assumption, by 
the function q(z). Determining q(z) from a comparison of the observed and Euclidean GRB burst rates 
one may determine therefore the z corresponding to a given C . For D = 4 the expression for z from 
q(z) is analytical, while for the remaining values of D it can be obtained numerically. The values of z 
corresponding to each C thus determined from the observations should, in a true unbounded standard 
candle flat model, be related to the photon emission rate C through 

C = 4irR 2 C[(l + z) 1/2 - l] 2 = 4 x 10 57 /i" 2 C*[(l + z) 1/2 - l] 2 photons s _1 , (12) 

where h = if o /(100 km s _1 Mpc -1 ), and where C is in units of photons/(cm 2 s). (If the average photon 
has an energy of e MeV, then L = 6.4 x 10 51 Ch~ 2 e[(l + z) 1 / 2 — l] 2 ergs/s is the peak bolometric luminosity 
in energy units). This procedure can be carried out for all C > Cth and for different D. If the results 
gave a constant C (or L ) for a given D, one would conclude that the BATSE data are compatible with 
the simplifying assumptions of the above model. 

The results of this exercise, when applied to the publicly available BATSE (2B) data are indicated in 
Table 1 for some specific values of C (c.f. also Figure 1). Values of C above 5 photons cm -2 s _1 were 

not included, since here the departure from C~ 3 I 2 is difficult to determine with the present approximate 
graphical eye-fit. From an inspection of the numbers in Table 1 it follows that there appears to be no 
drastic difference between the cases D = 3 and D ^ 3. A qualitatively similar conclusion has been 
reached, based on numerical integrations of the integrated burst counts, by Dermer (1992) and Piran 
(1992). One might conclude, very roughly, that C ~ const, for any D when C is in the interval 
(0.3-few) photons/(cm 2 s). Nevertheless, one gets systematically smaller values of C for C = (2 — 5) 
photons/(cm 2 s). For these large C the unsatisfactory fit of the theoretical curve of equ.(10) with the 
BATSE (2B) data may be caused by the assumption of a standard candle, it could be that the function 
q(z) is systematically larger than assumed here, or it may be that for C = (2 — 5) photons/(cm 2 s) a 
better treatment of the low number statistics and errors is needed. To decide between these options 
would require going beyond the present analytical and qualitative considerations, e.g. a more detailed 
statistical fit. This is currently in preparation, and will be presented elsewhere. However, it would not 
be surprising if the basic premises of the model in this Section are too simple. In the next Section we 
consider a more complete model, by relaxing the assumption of a standard candle assumption. 



4. THE EFFECT OF A LUMINOSITY FUNCTION 



Consider now a specific example of a photon emissivity distribution similar to a Schechter function, 

W = nC- m ) n {C/C min )-P , (13) 

where £ m8n < C < C max , and outside of this interval $(£) = 0. We assume that n = const and [3 is a 
dimensionless real constant. For the rest of this section we will further assume that fix) = X; i- e - we 
take a flat universe. 

Consider a given peak count C > Cth- Then the total number of GRB seen per year that are nearer 
than the comoving distance xi an d have peak emissivities between C and (C + dC) is given by (see equs.(5, 
7-9)) 

N'(< X i,C) = 47r$(£)rf£ r Rl(l - xf (A - D) X 2 d X 

Jo 
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4 8 ~ 2D 

= ^Rmc)dc x \ a *xj], (14) 
6 k=i 

where is given by (9). These bursters will have count rates larger than the value C corresponding to 
the emissivity C at distance xi, related through \i = (1 + Roi^C / C) 1 ! 2 ) -1 . Repeating this procedure 
for any £ m8n < C < C max for the same C , i.e. keeping C fixed but varying C and hence xi (and assuming 
a brightness limited sample; for a sample with a finite z max see Appendix B) one obtains 



where 



K x 3 ' 2 -?dx r (bx 1 l 2 f 
[1+ V a*- V 7 



! (l+&X 1 /2)3 L ^ K (1 + 6^1/2)^-1 

b*-*>J b (TT^ [1+ ^ a& (TT^F J; (16) 

1/2 

£ - b* 1 ' 2 , A=^»l; &= _Ah2_. (17) 



The integral 7 is a dimensionless non-negative number, depending on b (i.e. on C), D and on K. Then 
N(> C), as the function of C , depends on parameters n, C m i n , C max , [3 and D. 

If < b < A" 1 / 2 < 1, i.e. C> C e = £ max /(4TrR 2 ), 

I c f x^dx c ( \ K *' 2 - P ~ - V> 5 / 2 : (18) 

J i I m A, for j3 = 5/2 , 

since here one has always y < 1 with t/ <C 1 over most of the range; hence one can approximately take 
1 + y ~ 1 and one can also take the square bracket in (16) to be unity. Then from equ.(15) it follows that 

n ., 4tt C^ n n f (f - P)-^[{C max /C min )^ - 1] oc C- 3 ' 2 , for/3^f; 
l>OJ - 3 (4ttC7)3/2 \ ln(£ mM /£ mi „) oc C" 3 / 2 , for/3=§, ^ 

which is the Euclidean limit. 

If A' -1 / 2 < b < 1, or C/ = £ mjn /(47ri? 2 ) < C < C max / [Ait R 2 ) = C e , one may approximately write 

I „ 2b~ 5 + 2 ? $1 y 4 - 2 ?dy + 2A D b- 5 + 2 ? ^ y^dy 

{(2/(2/3 -5)) = const., for /? > 5/2; 

- In & 2 + (2A c /3) ~ const., for /? = 5/2; 

&- 5 + 2 / 3 (2/(5-2/3) + A c /(/3-l)), forl</3<5/2; (20) 

&- 3 ((2/3) + A c ln(& 2 A)), for/3=l; 

&- 3 A c A 1 -' 3 (l - for/3<l, 
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where 



Ad= [l + E2:? D a*]=3/ 1 (l-x) 8 - 2 Vdx 

= [3/(9 - 2D) - 6/(10 - 2D) + 3/(11 - 2D)] , (21) 

with cik given in (9). (Some particular values of Ad are A4 = 1, A3 = 1/10, A2 = 1/35, etc. For most 
of the range y < 1 one may use (1 + y) ~ 1, while for y > 1 one can take (1 + t/) ~ t/. For /? > (5/2) the 
second integral is much smaller than the first one, because Ad < 1 and b < 1. For [3 < 1 the first integral 
is much smaller than the second one, because b 2 K > 1. Roughly one may also take In b ~ const. In fact, 
one can check the validity of this approximation for several [3 by direct integration. For example, for 
13 = (3/2) and D = 3 one has exactly for b < 1 and &A' 1 / 2 > 1 7 = 0.72&" 2 . Equ.(20) gives I = l.2b~ 2 .) 
Substituting the value of I from equ.(20) into equ.(15) one obtains 

N(> C) ~ (47rn/3) x 

f ^n(/? - 5/2)- 1 (4ttC7)- 3 / 2 <x C" 3 / 2 , for /? > 5/2; 

^SL(^C)- s l 2 [(2A D li) - HC min /{R 2 ATtC))] oc C7- 3 / 2 , for /? = 5/2; 

^7„ 1 ^- 2,3 (4'rC') 1 -' 3 [(5/2-^)- 1 + (A i ,/(^-l))]ocC' 1 -' 3 , for 1< /3 < 5/2; (22) 

i? 3 [A c ln(£ ma:c /(47rC7/i 2 )) + (2/3)] ~ const, for /3 = 1; 

i? 3 ^(l - l3)- 1 {C max IC min ) 1 -l ) = const, for j3 < 1. 

A similar behavior for the range 1 < /? < 5/2 had been previously derived for a specific value of D = 3 
(Wasserman 1992), but here we see from equ.(22) that the scaling N(> C) oc C 1_l3 is independent of D, 
i.e. it is valid for an arbitrary redshift dependence of the comoving density. 

For ft > 1, or C< C/ = £ mjn /(47ri? 2 ) one may write 



,-K 

I ~ Ad b~ 3 / x'^dx, (23) 



where Ad is given in eq.(21). Hence 



N(> C) ~ (47r/3)ni? 3 ^ x Y ^ v~™*/~< 



f(P- 1)-![1 " {Cmax/Cmin) 1 -! 3 ] = Const. for {3 ± 1 ; 

for (3=1. 



(24) 



Equs.(19), (22) and (24) define the integral GRB count rate function N(> C) of the independent 
variable C; this function depends on the parameters h, D, C max , C m i n and [3. 



5. DISCUSSION 



The purpose of this paper is to provide analytical expressions characterizing the various regimes that 
would be present in a cosmological distribution of bursting sources with some generic luminosity and 
density distributions. With a simple redefinition of the density parameter D in eq.(2) to D ss = D + 1, 
the expressions are also valid for steady (non-bursting) sources. The specific example discussed is that 
of a cosmological distribution of gamma-ray burst sources. We do not present detailed statistical fits 
here, concentrating instead on extracting insights and quantitative estimates from the analytical models 
developed above, and analyze the effect of the physical assumptions and the parameters on the various 
asymptotic behaviors expected for the brightness distribution (log N- log C) of the bursts. 
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The simplest standard candle model is given here a general analytical representation, which comple- 
ments previous analyses by providing a specific dependence of the counts on the density evolution with 
redshift, as well as on the cosmological curvature (see App. A). This expression quantifies specifically 
the departures from a simple homogeneous Euclidean model as a function of either redshift or count 
rate. The log#(> C) — logC curve has two asymptotic behaviors, the Euclidean C~ 3 I 2 for high C 
(nearby sources) and C° for very low C (distant sources) and D < 9/2. The flattening, caused by 
cosmological redshift effects, is very gradual, being given by the function q(z) (the square bracket in 
equ.(ll)). Approximate fits to the BATSE and PVO data are possible (c.f. also Mao and Paczyhski, 
1992, Dermer, 1992, Piran, 1992, Fenimore & Bloom, 1995), which are not strongly constrained by the 
data. The asymptotic value of N(> C) is given by equ.(10), and is a function of D, the parameter 
characterizing the density dependence of sources via eq.(2). This assumes brightness-limited sources, or 

Xi(Ctft) = (1 + [^nRlCth/ Co] 11 " 2 )' 1 < Xmax = (1 - (1 + Zmm)" 1 ' 2 ), where z pa? is the redshift at which 
the most distant (but unobservable) sources are. Otherwise, for a redshift-limited sample, the N(> C) 
also becomes constant at the C below those corresponding to the weakest sources at z max , but the lim- 
iting value of N(> C) is smaller than (10) (see Appendix B). One advantage of the present analytical 
treatment is that, given an ideal set of data and assuming the brightness-limited standard candle case of 
§3, a simple and direct estimate of a lower limit to z max could be obtained from eq.(ll) as a function of 
the density evolution exponents D. One would do this by finding the z values at which flattening starts 
to set in (see also Table 1). In practice, however (see Fig. 1), this flattening regime is not clearly seen at 
the lowest C considered, and the data are not "ideal". While the standard candle model gives reasonable 
numbers (L ~ 6 x 10 51 erg/s, z max *J 1, depending on D), the present approximate fit is not unique 
(see §3). A detailed statistical fit is needed in order to get more quantitative constraints on such models 
(e.g. Emslie & Horack 1994, Cohen & Piran, 1995). However, a number of instrumental uncertainties 
and incompleteness consideration must be taken into account, especially at low C, and the constraints 
so far are weak. 

In the presence of a power-law luminosity function (13), equations (18-24) show that the relation 
log#(> C) vs. logC can have three different asymptotic behaviors, instead of two. A simple physical 
interpretation of this behavior for homogeneous Euclidean space is discussed in Wasserman (1992), which 
also indicates what can be expected in cosmology. Here, we have obtained the analytic behavior for 
cosmologies with a density evolution that can depart from a constant comoving density (D ^ 3 as well as 
D = 3), and have evaluated these affects also for Sl < 1. The values of the constants in the expressions 
for the counts (eqs. (19), (22), (24)) depend on D for a brightness limited sample. For a redshift limited 
sample the constants are somewhat different and can depend explicitly on z max (e.g. Appendix B), and 
for non-flat cosmologies they can also depend on Sl (see App. A). The slopes of (log#(> C)) increase 
from —3/2 at high C towards at low C , as in the standard candle case, but now for a luminosity function 
index 1 < [3 < (5/2) the integral counts go through an intermediate linear behavior where the slope is 
(1 — /?), before flattening to (for D < 9/2) at very low C, i.e., the behavior can be approximated by 
three straight lines rather than two. 

The first transition, from the Euclidean C~ 3 I 2 behavior to the C 1- ' 3 behavior, occurs at 

C * £gg.. (2.5, 

where Q e ~ [l + ^((5/2) - /?)/(/? - l)] 1/(/3_(5/2)) ~ l. (For /? > (5/2) the behavior C~ 3 ' 2 continues 
until the next transition Cj below; thus here no analogy of C e exists. For [3 < 1 equation (25) may 
also be used to define C e with Q e ~ (Ad((5/2) — — /3) -1 ) -2 / 3 ; here C e defines the break between 
the C~ 3 I 2 and C° behavior.). The observed log N — log C curve (Fig. 1) can be used for deriving an 
approximate estimate of the first break, C e ~ 6.9 photons cm -2 s _1 , which from eq.(25) yields 

C max ~ 2.5 x 10 58 /r 2 (C e /6.9) s" 1 . (26) 
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This corresponds, assuming a typical photon energy e MeV, to a luminosity L max ~ 2.5xl0 52 e/i 2 (C e /6.9) 
erg s _1 (c.f. also Wasserman, 1992). 

The second transition from the C 1_l3 behavior to the flat asymptotic behavior C° occurs at 

c ' = tl«" <*> 

where Q f ~ [l + (/? - l)((5/2) - /J) -1 ^ 1 ] 1/(/3_1) for 1 < /3 < 5/2. (Equ.(27) may be defined also for 

(3 > 5/2, with Qj ~ [(/? - l)/([/3 - 5/2]A c )] 2/3 . In this case Cj defines the break between the C*" 3 / 2 
and C° behavior. No analogy of C j exists for j3 < 1.) This second break point cannot be clearly defined 
from the data in Fig. 1, but an upper limit to it may be taken to be Cj & 0.4 photons cm -2 s _1 . (We 
note that the 2B catalogue has some data points at C & 0.4 photons cm -2 s _1 , for which however the 
triggering corrections become very uncertain, so we did not use any data points below this value). 

The specific value of [3 is most important in the range 1 < [3 < 5/2 because in this case there can be a 
range of C intermediate between the regimes C~ 3 I 2 and C° where the behavior of N(> C) reflects directly 
the effects of the luminosity function sampled at the largest redshifts. (Note that even for z max — ► oo 
as in §4 the relevant largest redshifts are z ~ 3 because that is where C e ~ £ ma:c /(47ri? 2 ) is equivalent 
to C(C max ) = Cmax/{^dl(z max )) = £ max /(4irR 2 [ v / l + z- l] 2 )). On the other hand, for j3 > 5/2 the 
N(> C) curve is dominated by the faintest bursts, while for (3 > 1 it is dominated by the brightest 
bursts in the luminosity function. From the data in Fig. 1, one can obtain an approximate fit to the 
intermediate slope of (1 — /?) ~ -0.88, or /3 ~ 1.88 (c.f. also Meegan et al. 1992, Wasserman, 1992). 
Note, however, that the credible contours for parameters of a wide variety of models are so far extremely 
large (e.g. Emslie & Horack, 1994, Cohen fePiran, 1995), so there is a degeneracy of possible models. 
For instance, Loredo and Wasserman (1992) make the point, also apparent from our analysis, that it is 
difficult to distinguish between a standard candle with no evolution and ones with a luminosity function, 
or standard candles with density evolution. A recent discussion of some of the relevant issues of data 
analysis methodology and incompleteness problems has been given by Loredo and Wasserman (1995). 

If one uses a fit with a power-law slope of 1 — [3 ~ —0.88, and assuming for the sake of argument D = 3 
so that A3 = 10 _1 and Qf ~ 2.2 x 10 1 , an upper limit for £ m8n is 

Cmin & 7 x W'-'h- 2 (C s /{)A)(Q s /2'i)- 1 s-\ (28) 

corresponding to a luminosity L m i„ ^ 7 x 10 49 £/i -2 (C//0.4)(Qj/22) -1 er g s_1 - The corresponding lower 
limit on the ratio of the maximum to minimum photon emissivity in the distribution function (13) is 

K = (Cmax/Cmin) * 3 x 10 2 (C e /6 .9) (C f /0 .4) " 1 (Q f /22) . (29) 

From Figure 1 we also obtain a lower limit value of N(> C j) ^ 3 x 10 2 averaged over approximately two 
years which, from eq.(23) with D = 3, Ad = 0.1 yields a lower limit to the total burst rate density 

n £ 1.5 x (A D /0A)- 1 h 3 Gpc- 3 yr- 1 . (30) 

If we assume galaxies to have a density of 0.01 h 3 Mpc -3 , this gives a lower limit to the total 
burst rate of TZ ^ 1.5 x 10~ 7 (A D /OA)' 1 bursts per year per galaxy. (Note that, from eqs.(3) and 
(13), n ~ n ). This lower limit can be checked also using equ.(19) at C e , which gives the same limit 
TZ <^ 1.5 x 10 _7 (A'/300)' 3_1 per galaxy per year. The K 13-1 dependence arises because of the assumed 
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C 1_l3 dependence with C oc C in this regime. The fact that these two independent estimates are similar, 
based on (23) evaluated for Ad = 0.1 and (19) which is independent of Ad, suggests that D = 3 is not 
an inconsistent assumption. 

The same analysis of §3-4 can be done for Sl ^ 1 (see Appendix A). The result is that small higher order 
terms appear in the integrals (8) and (14) which in (11) and (16) lead to additional small, Sl dependent 
correction terms in the square bracket, in addition to the previous small evolution and luminosity function 
correction terms, the leading term still being unity as before. These extra terms explicitly quantify the 
weakness of the dependence on Sl , and explain why numerical statistical fits (e.g. Loredo & Wasserman, 
1992, Emslie & Horack, 1995, Cohen & Piran, 1995) have not so far been able to constrain cosmological 
parameters. The reason for this insensitivity to the closure parameter is that the density, luminosity 
and cosmological effects are first order corrections in the comoving maximum distance xi (with Xi <C 1 
(App. A) for any phenomena dependent on star formation at z max & few). For different cosmological 
models the shape of the curves remains similar, and only the absolute length scales change, i.e. the way 
X depends on z. The latter could affect the observationally determined value of C max , C m i n (which also 
depend on h = (if o /100), the value of the present Hubble constant in units of 100 Km/s/Mpc, and thus 
on Sl ), but it does not allow one to discriminate between distance effects and possible intrinsic physical 
effects on C at the source. 

Another simplification made in §3 and §4 is to take the sources to be brightness-limited, with the 
density dependence n(z) oc (1 + z) D out to the largest z or smallest C observed. Realistically, however, 
one could expect a different functional dependence, or alternatively it might be that n(z) oc (l + z) D only 
for z < z max , while for z > z max one might have n(z) ~ 0. The straightforward generalization of the 
mathematical formalism of this paper to such redshift-limited cases is carried out in Appendix B. This 
shows that, while for z max & 3 the values of n would be somewhat affected, most effects associated with 
a finite z max are generally small. 

The luminosities and burst rate densities derived from this analysis are compatible with those expected 
for cosmological models based either on compact stellar merger or collapse scenarios (e.g. Narayan, 
Paczyhski & Piran 1992, Woosley 1993). Definitive conclusions regarding this would require further 
detailed statistical fits, and the resolution of a number of uncertainties about the completeness and the 
exposure function at low fluences. Models with a broad range of Sl , density evolution and luminosity 
function parameters appear acceptable. In particular, the present evidence does not appear sufficient to 
discriminate between a standard candle distribution (which would imply a well defined density) and a 
power law (or approximately similar) luminosity distribution (where the total density would be harder 
to estimate due to fading of counts at the low fluence end). However, some constraints can be made 
for specific cases. Thus, if the low fluence end of the counts is determined by a luminosity function 
of index [3 ~ 1.88 one has approximately 10% fewer sources in each decade of luminosity than in the 
previous lower decade, i.e. if K ~ 10 2 one has 90% of the sources between £ m8n and 10£ m8n and 10% 
between 10£ m8n and 10 2 £ m8n ~ Umax, a statement compatible with the findings of Emslie & Horack 
1994 and Ulmer & Wijers 1994. We note, however, that £ m8n (eq.[28]) is a current upper limit, which 
could become lower with more sensitive observations. While this cannot be ruled out energetically, if the 
integrated burst rates continued with a slope 1 — [3 ~ —0.88 down to count values much less than one 
order of magnitude below the present estimated upper limit Cf ~ 0.4 cm _2 s _1 , the required burst rate 
would exceed the maximum value 1Z ~ 10 -6 — 10 -5 per year per galaxy expected from compact mergers 
or collapse, providing constraints on such scenarios. 
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APPENDIX A: EFFECTS OF A NON-FLAT COSMOLOGY 



For a general Friedmann model the integrated number of bursts out to xi is still given by eq.(5), but 
f(x) is sinhx (sinx) for an open (closed) universe. The other difference is that the expansion factor 
(which was R(r])/R = r/ 2 for the flat universe) is now, assuming Sl < 1, 



R(r)) (cosh r) — 1) 1 



(31) 



R (cosh r] — 1) 1 + z 
where cosh rj = (2/Sl ) — 1, x = Vo ~ V an( i cosh rj = cosh(ry — x)- F° r small Sl we can expand 

cosh r] - 1 ^ l _ 
cosh T] — 1 ' 

where g ~ [1 + -|S! ] ~ 1. Also, since for z\ ^ few one has xi "C 1, we can take 

/ 2 ( X ) = smh 2 x~X 2 . (33) 
With this we have the integrated number of bursts per year out to xi as 

N'(< xi) ~ ^(C)dCRl Jf (1 - gx) 2(4 - D) X 2 dx 

~ f R 3 M£)djC xl [1 + El ( ^ D) a k (g X i) k ] , (34) 

where aj, is given by (9). Since Xi <C 1 and g ~ 1, the square bracket is again [• • •] ~ 1 and we get the 
same result as eq.(14) (to within small correction terms of order g -1 ~ 1 — (Sl /2) — ► 1 for Sl ^ 1). 
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Therefore N(< C) is also given by eq.(16) to within the same small correction factor, and (19), (21) 
and (24) are also approximately correct for Sl <C 1. The specific dependence of x on z is, however, a 
function of Sl , which affects the length and luminosity scales, but not the shapes of the curves. Other 
effects associated with a possible cosmological constant can also be included in a straightforward manner 
(e.g. Piran 1992, Cohen, Kolatt & Piran 1994, Emslie & Horack 1994), and have been ignored here for 
simplicity. 



APPENDIX B: THE EFFECT OF A LIMITING REDSHIFT 



Consider the case when the density dependence of equ.(3) is valid only up to some maximum redshift 



z. 



max , 



n(z) = n (l + z) D (l-e(z/z max )) , (35) 
n = f£-*(C)dC = nJ 1 K x-<>dx= J (1 " W^) 1 "' 3 - 1], for ^ 1, 

Z """> 1 [l*l(£ma X /£min), for /? = 1 , 

where 0(z / ' z max ) is the Heaviside function. There is in this case a limiting comoving radial coordinate 
Xmax defined unambiguously by z max . (For the flat universe f(x) = X one has Xmax = (l-[l + 2maj~ 1/2 ), 
and in what follows we restrict ourselves to this case only.) Equ. (5) is correct only for xi < Xmax- For 
Xi > Xmax one has to take as upper limit of integration in equ. (5) Xmax instead of xi- Equs.(6-7) are 
again correct, but, of course, are defined only for z < z max . 

The standard candle model considerations of §3 can be generalized as follows. Equs.(8-9) are fulfilled 
only in the case when xi < Xmax, corresponding to C > C'f- The value C'f is defined by Xmax = 
(1 + R [At{Cf I C-o] 1 ^ 2 ) -1 ■ For xi > Xmax, in equ.(8) the upper limit of integration is Xmax- This 
means that for C > C'f, equ. (9) describes as before the function N(> C), but for C < C'f one has 
N(> C) = N(> C'f) = const, (there is a horizontal straight line for C < C'f)- This also means that in 
equ.(10) the appropriate limit for C — ► is smaller than the value given, and in equ.(ll) the function 
q(z) is also defined only for z < z max . 

The generalization of the power-law luminosity function model of §4 may be done as follows. Here 

again equ. (14) holds only for xi < Xmax- When xi > Xmax = (1 - (1 + Zmaj;)" 1/2 ), the upper limit 
of integration is Xmax- For an arbitrary C in the interval < C < oo one can define a peak photon 
emissivity C(C') and a critical count rate C'e through 



C(C)= C47rR 2 oX 2 max(l ~ Xmax)- 2 = GAitR 2 ^! + z max - I f = C47rdl(z max ), 
C E = £max/(4:Trdl(z max )) , (36) 



where d c (z max ) = R (\/l + z max — 1) is the count distance of eq. (7), Ce is the finite z equivalent of C e 
and < jC(C) < oo. 

Consider first the case when C(C) > C max , i.e. C > Ce (or b < [i/l + z max — 1]A' -1 / 2 ). This 
means that for any C, where £ m8n < C < C max < C{C), the corresponding comoving coordinate xi 
defined by xi = (1 + R {^C / C) 1 ! 2 ) -1 is smaller than Xmax, and therefore in this case the considerations 
of §4 are correct without any change. In other words, for C > Ce the function N(> C) vs. C is 
unchanged. (Some numerical examples are useful here. For example, for Xmax = 0.9 (z max = 99), one 
has Ce = (l/81)£ ma j;/(47ri? 2 ). Then, if A' < 81, equs. (18-21) remain unchanged. But, if Xmax = 0.5 
(z max = 3), one has Ce = £-max I '(4tt 'R 2 ) = C e ss defined in equ. (25), and only the Euclidean part of 
equs. (18-19) remain unchanged. For z max < 3 even the Euclidean part is influenced by the existence of 
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Zmax- Nevertheless, the regimes C 3 I 2 , C 1 @ and C° still occur, but in modified form and with different 
normalizations and break points; see below). 

Next, take C min < C{C) < C ma x, or C F = £ m i„/(4:ird 2 c (z max )) < C < C max / {^d 2 .(z max )) = C E 
(or, equivalently, [VI + z max ~ l]K~ 1/2 < b < [VI + z max ~ 1])- For C min < C < C(C) equ.(14) still 
holds, because here xi < Xmax- On the other hand, for C{C) < C < C max in equ.(14) the upper limit of 
integration is Xmax- Hence, using the distribution (13), it follows that 

NO C) - 4 - £ ^ r )IC - X " Pdx [1 + 8 V a {bxl ' 2)k 1 
1 ' ~ 3(4ttC7) 3 / 2 J, (1 + bx 1 /*) 3 [ fr[ (1 + bx^ f J 

4tt 8 ~ 2D f K 

A k=1 Jc(c)/c m , n 

= yi? 3 n b 3 h + yi? 3 n B(z max ,D) I 2 , (37) 

where b is defined by equ.(17), C{C) is defined by equ.(36), and B(z max ,D) is a constant. We see that 
N(> C) is defined by the sum of two terms having different behaviors. The first term is practically 
identical to the relation defined by equs. (15-16); the only difference follows from the upper limit of 
integration in I\. The second term is exactly calculable; the dependence on C follows from the fact that 
the lower limit of integration in I 2 depends on C via C{C). The integral I2 is given by 

K x ~P dx= {(l-P)-%Crnax/C mln Y- p -(C(C)/C mln Y- p l for /? ^ 1 , ( g g) 

k \\n(jC max /jC(Cj), for /? = 1, 

where, for brevity, we have denoted K = (£(C)/£ m8n ). The integral I\ may be calculated as in §4. 
There are three different cases. First, if z max < 3 (or Xmax < 0.5), then b < 1, and, in addition, 

bx 1 ! 2 = (xmax/(l - Xmax)) = (VI + z max ~ 1) < 1 for a; = C(C)/C min . In this case one may write 

{(b/2-p)- l [C(C)/C min fl 2 -^ -I] for fit 5/2, . , 

1 " \ln(£(C)/£ mi „), for = 5/2. K ™> 

Second, consider the range 3 < z max < ((1 + A' 1/2 ) 2 - 1) (or 0.5 < Xmax < A' 1/2 /(l + A 1/2 )). In this 

case both b < 1 and b > 1 may occur; but always bx 1 ! 2 = Xmax / '(1 — Xmax) > 1 for x = C(C)/C m i n . 
Therefore, if b < 1, one has 



h 



(2/(2/3- 5)) = const., for /? > 5/2; 

- In b 2 + (2A D /3) ~ const., for /? = 5/2; 

&- 5 + 2 / 3 (2/(5-2/3)+A c /(/3- 1)), forl</3<5/2; (40) 

&- 3 ((2/3) + ^ln(& 2 A')), for/3=l; 

6" 3 yl z ,A' 1 -^(l - for/3<l, 
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where we have proceeded as in the derivation of (20). The results are also the same, except for the fact 

that K is substituted by K for [3 < 1. In the same range of z max but for b > 1, instead of eq.(40) one 
obtains 

h - A D h-* f .6-3 x { fMl - (^.-n) 1 -"], for £ # 1, (41) 

J i L ln(£/£ mi „), for (i = 1. 

Third, consider z ma:c > ((1 + A' 1 / 2 ) 2 — 1) (or Xmai > A' 1 / 2 /(l + A' 1 / 2 )). In this case one necessarily has 
& > 1, and 7i is again given by equ.(41). Substituting I\ and I2 into equ.(37) one obtains the function 
N(> C) vs. C, which depends now on the parameters £ m in, C max , [3, h, z max , and D in analytical 
form within this range of C . 

Finally, consider the case when C(C) < jC m i„; i-e. C < Cf = £-minl{^'xd 2 c (z, max )), with c? c given by 
eq.(7), and Cf is the finite z equivalent of Cf. In this case, one always has xi > Xmai- Therefore, here 
equ.(14) does not hold, because the upper limit of integration is Xmax- Hence one obtains 

47T 8-2£> ~k ^ 

N(> C) = -^-R 3 n Xmax [1 + J2 ak x ^ ] / x ~ Pdx = -^ R l n oB{z max , D) = const, (42) 
6 k=i Jl 6 

where n is given by equ.(35) and B(z max ,D) by equ.(37). This means that here, as in the standard 
candle case, there is again a horizontal straight line on the graph N(> C) vs. C . (The differential counts 
N(C) at C < C'f are expected to be zero, while the integral counts N(> C) of (42) are constant and 
given by the value at C = Cf)- 

These relations would allow one in principle to determine the parameters C m i n , C max , [3, D, n and 
Zmax- For instance, [3 could be determined similarly to the case with z max = 00 with large C from the 
asymptote of the slope (if 1 < [3 < (5/2)). Then from the transition in the graph N(> C) vs. C between 
the C~ 3 I 2 to the C 1- ' 3 behavior one can determine C max (if z max > 3 equ.(26) holds). Using this, from 
the Euclidean part, one obtains n (see equs.(15) and (18)). In the ideal case then Xmax (or z max ) could 
be obtained from a careful analysis of the departure of the N(> C) graph from a straight line with slope 
(1 — f3), because this departure should be given by the second term in equ.(37), which depends on Xmax- 
Then the horizontal C° straight line at very small C , and its break to the C 1_l3 behavior, could yield the 
value of K (hence of £ m in) and D via equs.(37), (38) and (42). 

In practice, there would be severe difficulties due to incompleteness and errors at the lowest count 
rates. It may be very hard, if not impossible, to determine the C° straight line expected at the lowest C . 
Even in such a less than ideal N(> C) graph, the formulas of this Appendix would provide a strategy 
for obtaining bounds on z max . To illustrate this, consider the case when 1 < [3 < 5/2, and z max < 3. In 
this case the function N(> C) is given by 

N(> C) ~ (4tt/3)x 
(llcr/i g ~ 3/2 > for C>C E , 

r nC^L U(C)IC m , n Y'i-P-l] 

L(4ttC) 3 / 2 5/2-/3 (43) 

+^n5(z ma ,,J) ^-/ £ -^ 1 -;_-f^/ £ -^^ ], forCV <C<C E , 
B? n B(z max ,D) = const, for C < C F , 

where C E = C max / {^dl(z max )) , C F = Cmin / {^dl(z max )) are the equivalents of C e , Cf in §4, d c (z max ) = 
Ro(V^ + z max — 1), £{C) is defined by equ.(36), and n = h(l — /3) _1 [A' 1- ' 3 — 1]. We see that the Eu- 
clidean part is again continued by an approximate straight line with slope (1 — j3) for Cf ^ C & C E . 
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However, due to the second term in the square bracket of equ.(42), there is a departure from this straight 
slope (1 — /?) as C decreases. If the data allowed it, a careful analysis might determine this second term 
from this departure, and hence also z max . If, in addition, a C° portion were unequivocally discernible 
(which may not be the case), an additional relationship would be obtained linking z max , D and n . 

Table 1 
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0.08 
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0.23 


0.27 


0.34 
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0.45 


0.50 


0.55 


0.66 


C h 2 /W 56 


0.3 


0.4 


0.8 


1.0 


1.0 


1.0 


1.0 


1.0 


1.0 


1.0 
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Table Caption 



Table. 1 Comparison of the cosmological unbounded standard candle assumption with comoving densities 
depending as (1 + z) D for selected values of the count rate C . The case D = 3 describes the constant 
comoving density; D = 4 gives the case when the comoving density is increasing with z; D < 3 gives the 
cases when the comoving density is decreasing with z. If correct, one would obtain, for a given D, the 
same C at all z or C. Here C is in units of photons/(cm 2 s), C is in units of photons/s. 



Figure Caption 

Fig.l. Graph of log#(> C) versus C as obtained from the publicly available BATSE 2B catalogue. Two 
straight lines with slopes -3/2 and -0.88 are also shown superimposed, representing the Euclidean limit 
and a particular power-law luminosity function distribution (see text). 
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